Abstract
Background We have run the simplified Naomi model using a range of inference methods.
Task In this report, we compare the accuracy of the posterior distributions obtained from these inference methods.
tmb <- readRDS("depends/tmb.rds")
aghq <- readRDS("depends/aghq.rds")
tmbstan <- readRDS("depends/tmbstan.rds")
All of the possible parameter names are as follows:
unique(names(tmb$fit$obj$env$par))
## [1] "beta_rho" "beta_alpha" "beta_lambda" "beta_anc_rho" "beta_anc_alpha" "logit_phi_rho_x" "log_sigma_rho_x" "logit_phi_rho_xs" "log_sigma_rho_xs" "logit_phi_rho_a"
## [11] "log_sigma_rho_a" "logit_phi_rho_as" "log_sigma_rho_as" "log_sigma_rho_xa" "u_rho_x" "us_rho_x" "u_rho_xs" "us_rho_xs" "u_rho_a" "u_rho_as"
## [21] "logit_phi_alpha_x" "log_sigma_alpha_x" "logit_phi_alpha_xs" "log_sigma_alpha_xs" "logit_phi_alpha_a" "log_sigma_alpha_a" "logit_phi_alpha_as" "log_sigma_alpha_as" "log_sigma_alpha_xa" "u_alpha_x"
## [31] "us_alpha_x" "u_alpha_xs" "us_alpha_xs" "u_alpha_a" "u_alpha_as" "u_alpha_xa" "OmegaT_raw" "log_betaT" "log_sigma_lambda_x" "ui_lambda_x"
## [41] "log_sigma_ancrho_x" "log_sigma_ancalpha_x" "ui_anc_rho_x" "ui_anc_alpha_x" "log_sigma_or_gamma" "log_or_gamma"
data.frame(
"TMB" = tmb$time,
"aghq" = aghq$time,
"tmbstan" = tmbstan$time
)
histogram_plot("beta_anc_rho")
ks_helper <- function(par) to_ks_df(par) %>% ks_plot(par = par)
ks_helper("beta")
ks_helper("logit")
ks_helper("log_sigma")
ks_helper("u_rho_x")
ks_helper("u_rho_xs")
ks_helper("us_rho_x")
ks_helper("us_rho_xs")
ks_helper("u_rho_a")
ks_helper("u_rho_as")
ks_helper("u_alpha_x")
ks_helper("u_alpha_xs")
ks_helper("us_alpha_x")
ks_helper("us_alpha_xs")
ks_helper("u_alpha_a")
ks_helper("u_alpha_as")
ks_helper("u_alpha_xa")
ks_helper("ui_anc_rho_x")
ks_helper("ui_anc_alpha_x")
ks_helper("log_or_gamma")
ks_summary <- lapply(unique(names(tmb$fit$obj$env$par)), function(x) {
to_ks_df(x) %>%
group_by(method) %>%
summarise(ks = mean(ks), par = x)
}) %>%
bind_rows() %>%
pivot_wider(names_from = "method", values_from = "ks") %>%
rename(
"Parameter" = "par",
"KS(aghq, tmbstan)" = "aghq",
"KS(TMB, tmbstan)" = "TMB",
)
ks_summary %>%
gt::gt() %>%
gt::fmt_number(
columns = starts_with("KS"),
decimals = 3
)
| Parameter | KS(aghq, tmbstan) | KS(TMB, tmbstan) |
|---|---|---|
| beta_rho | 0.182 | 0.188 |
| beta_alpha | 0.117 | 0.117 |
| beta_lambda | 0.057 | 0.058 |
| beta_anc_rho | 0.111 | 0.110 |
| beta_anc_alpha | 0.097 | 0.057 |
| logit_phi_rho_x | 0.391 | 0.626 |
| log_sigma_rho_x | 0.581 | 0.656 |
| logit_phi_rho_xs | 0.270 | 0.638 |
| log_sigma_rho_xs | 0.839 | 0.646 |
| logit_phi_rho_a | 0.828 | 0.528 |
| log_sigma_rho_a | 0.438 | 0.557 |
| logit_phi_rho_as | 0.824 | 0.514 |
| log_sigma_rho_as | 0.294 | 0.542 |
| log_sigma_rho_xa | 0.438 | 0.692 |
| u_rho_x | 0.116 | 0.114 |
| us_rho_x | 0.103 | 0.102 |
| u_rho_xs | 0.147 | 0.146 |
| us_rho_xs | 0.101 | 0.101 |
| u_rho_a | 0.179 | 0.180 |
| u_rho_as | 0.125 | 0.118 |
| logit_phi_alpha_x | 0.317 | 0.571 |
| log_sigma_alpha_x | 0.709 | 0.592 |
| logit_phi_alpha_xs | 0.279 | 0.587 |
| log_sigma_alpha_xs | 0.686 | 0.613 |
| logit_phi_alpha_a | 0.791 | 0.533 |
| log_sigma_alpha_a | 0.290 | 0.540 |
| logit_phi_alpha_as | 0.749 | 0.511 |
| log_sigma_alpha_as | 0.237 | 0.540 |
| log_sigma_alpha_xa | 0.706 | 0.653 |
| u_alpha_x | 0.091 | 0.086 |
| us_alpha_x | 0.091 | 0.091 |
| u_alpha_xs | 0.078 | 0.069 |
| us_alpha_xs | 0.116 | 0.114 |
| u_alpha_a | 0.126 | 0.124 |
| u_alpha_as | 0.086 | 0.080 |
| u_alpha_xa | 0.073 | 0.072 |
| OmegaT_raw | 0.192 | 0.548 |
| log_betaT | 0.174 | 0.679 |
| log_sigma_lambda_x | 0.556 | 0.770 |
| ui_lambda_x | 0.205 | 0.205 |
| log_sigma_ancrho_x | 0.714 | 0.524 |
| log_sigma_ancalpha_x | 0.706 | 0.616 |
| ui_anc_rho_x | 0.066 | 0.071 |
| ui_anc_alpha_x | 0.132 | 0.134 |
| log_sigma_or_gamma | 0.507 | 0.521 |
| log_or_gamma | 0.075 | 0.075 |
ggplot(ks_summary, aes(x = `KS(TMB, tmbstan)`, y = `KS(aghq, tmbstan)`)) +
geom_point(alpha = 0.5) +
xlim(0, 1) +
ylim(0, 1) +
geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
labs(x = "KS(aghq, tmbstan)", y = "KS(TMB, tmbstan)") +
theme_minimal()
ks_summary %>%
mutate(`KS difference` = `KS(TMB, tmbstan)` - `KS(aghq, tmbstan)`) %>%
ggplot(aes(x = `KS difference`)) +
geom_boxplot(width = 0.5) +
coord_flip() +
labs(x = "KS(TMB, tmbstan) - KS(aghq, tmbstan)") +
theme_minimal()
#' To write!
#' To write!